/*
 * EvoAlghorithm.cpp
 *
 *  Created on: 15.11.2012
 *      Author: shadius
 */

#include "EvoAlghorithm.h"
#include <boost/random.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/thread.hpp>

EvoAlghorithm::EvoAlghorithm(const shared_ptr<IFastDataApprox> &BGValues,
		const shared_ptr<IFastDataApprox> &IGValues,
		const shared_ptr<vector<double>> p_ECGValues, INPUT_PARAMS params) :
		p_BGValues(BGValues), p_IGValues(IGValues), p_population(
				new boost::numeric::ublas::matrix<double>(POPULATION_SIZE,
						PARAMETER_COUNT)), p_fitness(
				new vector<double>(POPULATION_SIZE)), p_ECGValues(p_ECGValues), params(
				params) {

	srand(time(NULL) + getpid());

}

void EvoAlghorithm::mutate(double& mutationRate,
		shared_ptr<boost::numeric::ublas::matrix<double> >& p_newPopulation,
		int& iInd) {
	//mutace
	for (int child = 0; child < 2; child++) {
		for (int i = 0; i < PARAMETER_COUNT; i++) {
			if (i == PI && params.isPset)
				continue;
			if (i == CI && params.isCset)
				continue;

			if (random() < mutationRate) {
				double randomValue = 0.2 * sqrt(-2 * log(random(0, 1)))
						* sin(2 * M_PI * random(0, 1));
				//double randomValue = getRandomDoubleUsingNormalDistribution(0, 0.2);
				randomValue += (*p_newPopulation)(iInd + child, i);
				if (checkRange(i, randomValue)) {
					(*p_newPopulation)(iInd + child, i) = randomValue;
				}
			}
		}
	}
}
bool EvoAlghorithm::checkRange(int index, double value) {
	switch (index) {
	case PI:
		return checkRangeP(value);
		break;
	case CGI:
		return checkRangeCG(value);
		break;
	case CI:
		return checkRangeC(value);
		break;
	case DTI:
		return checkRangeDT(value);
		break;
	default:
		return false;
		break;
	}

}
double EvoAlghorithm::getClosestBound(double min, double max, double num) {
	double minLim = fabs(num - min);
	double maxLim = fabs(num - max);

	return (minLim < maxLim ? (min + 0.01) : (max - 0.01));

}
void EvoAlghorithm::initializePopulation() {
	//inicializujeme populaci
	for (int i = 0; i < POPULATION_SIZE; i++) {

		(*p_population)(i, PI) = (
				params.isPset ?
						params.p :
						random(PARAM_P_DEFAULT_MIN, PARAM_P_DEFAULT_MAX));
		(*p_population)(i, CGI) = random(PARAM_CG_DEFAULT_MIN,
				PARAM_CG_DEFAULT_MAX);
		//cout << "CG: " << (*p_population)(i, CGI) << endl;
		(*p_population)(i, CI) = (
				params.isCset ?
						params.c :
						random(PARAM_C_DEFAULT_MIN, PARAM_C_DEFAULT_MAX));

		double k = 0;
		double dt = 0;
		findKandDT(&k, &dt);
		(*p_population)(i, DTI) = dt;
		(*p_population)(i, KI) = k;
	}
}

void EvoAlghorithm::fitnessParalel(int i, int offset, LTHREAD_RESULTS results,
		LPAPPROXPOINT bgAPPROXPOINTS, LPAPPROXPOINT igAPPROXPOINTS, int len) {

	double localBestFitness = 0;
	double localbestFitnessIdx = 0;
	double localMinFitness = numeric_limits<double>::infinity();

	for (int j = i; j < i + offset; j++) {
		(*p_fitness)[j] = fitness(j, bgAPPROXPOINTS, igAPPROXPOINTS, len);

//	//	cout << "nalezena fitness " << (*p_fitness)[i] << endl;
		if ((*p_fitness)[j] > localBestFitness) {
			localBestFitness = (*p_fitness)[j];
			localbestFitnessIdx = j;
		}
		if ((*p_fitness)[j] < localMinFitness) {
			localMinFitness = (*p_fitness)[j];
		}
	}
	results->bestFitness = localBestFitness;
	results->bestFitnessIdx = localbestFitnessIdx;
	results->minFitness = localMinFitness;

}

RESULT EvoAlghorithm::compute(LPAPPROXPOINT bgAPPROXPOINTS,
		LPAPPROXPOINT igAPPROXPOINTS, int len) {

	THREAD_RESULTS results[params.threadCount];
	boost::thread threads[params.threadCount];
	double mutationRate = MIN_MUTATION_RATE;

	cout << "vypočteme fitness vsech jedincu v populaci" << cout;

	//inicializujeme populaci
	initializePopulation();

	double previousFitness = 0.0f;
	for (int iGeneration = 0; iGeneration < GENERATIONS_NUM; iGeneration++) {
		double bestFitness = 0.0f;
		int bestFitnessIdx = 0;
		double minFitness = numeric_limits<double>::infinity();

		std::shared_ptr<boost::numeric::ublas::matrix<double>> p_newPopulation =
				make_shared<boost::numeric::ublas::matrix<double>>(
						POPULATION_SIZE, PARAMETER_COUNT);

		int offset = POPULATION_SIZE / params.threadCount;
		int i = 0;
		int threadIndex = 0;
		while (i < POPULATION_SIZE) {
			if (i + 2 * offset > POPULATION_SIZE) {
				offset = POPULATION_SIZE - i;
			}
			threads[threadIndex] = boost::thread(
					(boost::bind(&EvoAlghorithm::fitnessParalel, this, i,
							offset, &results[threadIndex], bgAPPROXPOINTS,
							igAPPROXPOINTS, len)));
			i += offset;
			threadIndex++;
		}
		for (int i = 0; i < params.threadCount; i++) {
			threads[i].join();
		}
		for (int i = 0; i < params.threadCount; i++) {
			if (results[i].bestFitness > bestFitness) {

				bestFitness = results[i].bestFitness;
				bestFitnessIdx = results[i].bestFitnessIdx;
			}
			if (results[i].minFitness < minFitness) {
				minFitness = results[i].minFitness;
			}
		}

		//vypočteme fitness vsech jedincu v populaci
//		for (unsigned i = 0; i < POPULATION_SIZE; i++) {
//
//			fitnessParalel(i, bgAPPROXPOINTS, igAPPROXPOINTS, len, bestFitness,
//					bestFitnessIdx, minFitness);
//
//		}

// upravime hodnoty fitness pro lepsi selekci od 0 do bestFitness-minFitness
		double sumFitness = 0.0;
		for (int i = 0; i < POPULATION_SIZE; i++) {
			(*p_fitness)[i] -= minFitness;
			sumFitness += (*p_fitness)[i];
		}

		// vytvoreni nove populace, ktera nahradi stavajici
		for (int iInd = 0; iInd < POPULATION_SIZE; iInd += 2) {
			int individuumIdx1 = pickRandomIndividuum(sumFitness);
			int individuumIdx2 = pickRandomIndividuum(sumFitness);
			if (individuumIdx1 < 0 || individuumIdx2 < 0) {
				cout << "to se mi nezdá" << endl;
			}
			//prirozeny vyber
			if (random() < CROSS_RATE) {
				for (int i = 0; i < PARAMETER_COUNT; i++) {
					if (random() < 0.5) {
						(*p_newPopulation)(iInd, i) = (*p_population)(
								individuumIdx1, i);
						(*p_newPopulation)(iInd + 1, i) = (*p_population)(
								individuumIdx2, i);

					} else {
						(*p_newPopulation)(iInd, i) = (*p_population)(
								individuumIdx2, i);
						(*p_newPopulation)(iInd + 1, i) = (*p_population)(
								individuumIdx1, i);
					}
				}
				//mutace
				mutate(mutationRate, p_newPopulation, iInd);

			} else { //neni krizeni
				for (int i = 0; i < PARAMETER_COUNT; i++) {
					(*p_newPopulation)(iInd, i) = (*p_population)(
							individuumIdx1, i);
					(*p_newPopulation)(iInd + 1, i) = (*p_population)(
							individuumIdx2, i);

				}
			}
		}

		//elitismus - zachování nejlepšího
		for (int i = 0; i < PARAMETER_COUNT; i++) {
			(*p_newPopulation)(0, i) = (*p_population)(bestFitnessIdx, i);
		}
		p_population.reset();

		p_population = p_newPopulation;
		if (bestFitness > previousFitness) {
			cout << "Generation " << iGeneration << ", best fitness "
					<< bestFitness << endl;
			previousFitness = bestFitness;
			mutationRate = MIN_MUTATION_RATE;
			//V první generaci začít s nejnižší mutací a v každé následující generaci mutaci zvýšit.
			//Mutace se znova vrátí na nejnižší hodnotu v případě, že v nové generaci je nalezen lepší jedinec
			//než v předchozí anebo v případě, že mutace dosáhla své maximální hodnoty.
			//
			//
		} else {
			if (mutationRate == MAX_MUTATION_RATE)
				mutationRate = MIN_MUTATION_RATE;
			else
				mutationRate += MUTATION_STEP;

		}
	}
	cout << "p = " << (*p_population)(0, 0) << endl;
	cout << "cg = " << (*p_population)(0, 1) << endl;
	cout << "c = " << (*p_population)(0, 2) << endl;
	cout << "k = " << (*p_population)(0, 3) << endl;
	cout << "dt = " << (*p_population)(0, 4) << endl;
	double sum = 0;
	int count = 0;
	double p = (*p_population)(0, PI);
	double cg = (*p_population)(0, CGI);
	double c = (*p_population)(0, CI);
	double dt = (*p_population)(0, DTI);
	double k = (*p_population)(0, KI);
	if (countAverage(bgAPPROXPOINTS, igAPPROXPOINTS, p, cg, c, dt, k, len, &sum,
			&count) == E_INVALIDARG) {
		cout << "Pruser!!" << endl;
		for (unsigned i = 0; i < p_ECGValues->size(); i++) {
			if (!checkTimeJumpRange(k, dt, (*p_ECGValues)[i])) {
				cout << "FAIIIIL JUMPY!!" << endl;
			}
		}

		if (!checkRangeDT(dt)) {
			cout << "FAIIIIL DT!!" << endl;
		}
	} else {
		cout << "soucet rozdilu je: " << sum;
		cout << "pocet prvku je: " << count;
		cout << "average je: " << sum / (double) (count);

	}

	RESULT res;
	res.c = c;
	res.cg = cg;
	res.dt = dt;
	res.k = k;
	res.p = p;

	return res;
}
int EvoAlghorithm::pickRandomIndividuum(double sumFitness) {
	double sum = 0.0;
	double fitnessPoint = random() * sumFitness;
	for (int i = 0; i < POPULATION_SIZE; i++) {
		sum += (*p_fitness)[i];
		if (sum > fitnessPoint)
			return i;
	}
	return 0;
}

double EvoAlghorithm::fitness(int iInd, LPAPPROXPOINT bgAPPROXPOINTS,
		LPAPPROXPOINT igAPPROXPOINTS, int len) {
	double p = (*p_population)(iInd, PI);
	double cg = (*p_population)(iInd, CGI);
	double c = (*p_population)(iInd, CI);
	double dt = (*p_population)(iInd, DTI);
	double k = (*p_population)(iInd, KI);

	int count = 0;
	double sum = 0;
	if (countAverage(bgAPPROXPOINTS, igAPPROXPOINTS, p, cg, c, dt, k, len, &sum,
	//TODO:něco udělat s tím E_INVALIDARG
			&count) == E_INVALIDARG) {
		return 0;
	}
	return 1.0 / (sum / (double) count);

}

HRESULT EvoAlghorithm::countAverage(LPAPPROXPOINT bgAPPROXPOINTS,
		LPAPPROXPOINT igAPPROXPOINTS, double p, double cg, double c, double dt,
		double k, int len, double *sum, int *count) {

	double localSum = 0;
	int localCount = 0;
	for (int i = 0; i < len; i++) {

		int igPoint = 0;
		double ecg = (*p_ECGValues)[i];
		if (!checkRangeC(c) || !checkRangeCG(cg) || !checkRangeP(p)
				|| !checkRangeDT(dt) || !checkTimeJumpRange(k, dt, ecg)) {
			return E_INVALIDARG;
		}
		p_BGValues->GetPrecedingPoint(
				TimeUtils::convertSecToWin(i) + dt + k * dt * ecg, &igPoint);

		if (igPoint < len && igPoint > 0) {
			double dif = fabs(
					p * bgAPPROXPOINTS[i].y + cg * bgAPPROXPOINTS[i].y * ecg + c
							- igAPPROXPOINTS[igPoint].y);

			localCount++;
			localSum += dif;
		}

	}
	*count = localCount;
	*sum = localSum;
	return S_OK;
}

bool EvoAlghorithm::checkRangeP(double p) {
	if (params.isPset) {
		return p == params.p;
	}
	return p > PARAM_P_DEFAULT_MIN && p < PARAM_P_DEFAULT_MAX;
}

bool EvoAlghorithm::checkRangeDT(double dt) {
	return dt > params.dtMin && dt < params.dtMax;
}

bool EvoAlghorithm::checkRangeCG(double cg) {
	return cg < 0;
}

bool EvoAlghorithm::checkRangeC(double c) {
	if (params.isCset) {
		return c == params.c;
	}
	return c > 0;
}

bool EvoAlghorithm::checkTimeJumpRange(double k, double dt, double ecg) {
	double jmp = k * dt * ecg;

	return jmp <= 5 * OneMinute && jmp >= -5 * OneMinute;
}

double EvoAlghorithm::random() {
	return ((double) rand() / (double) RAND_MAX);
}

double EvoAlghorithm::random(double min, double max) {
	double r = random();
	return min + r * (max - min);
}

//double EvoAlghorithm::randomGauss() {
//	std::mt19937 generator;
//	double mean = 0.0;
//	double std  = 1.0;
//	std::normal_distribution<double> normal(mean, std);
//	for(int i = 0; i < 50; i++)
//	cerr << "Normal: " << normal(generator) << endl;
//	return 0.0;

//}

void EvoAlghorithm::findKandDT(double *k, double *dt) {
	vector<double>::const_iterator it;

	double found = true;
	double localK;
	double localDt;
	while (found) {
		if (params.isKset) {
			localK = random(params.kMin, params.kMax);
		} else {
			localK = random();
		}
		localDt = random(params.dtMin, params.dtMax);

		for (unsigned i = 0; i < p_ECGValues->size(); i++) {
			found = !checkTimeJumpRange(localK, localDt, (*p_ECGValues)[i]);
			if (found)
				break;
		}
	}
	*k = localK;
	*dt = localDt;

}

EvoAlghorithm::~EvoAlghorithm() {
}
